Multistate Bennett Acceptance Ratio Method Research Articles

Bayesian Multistate Bennett Acceptance Ratio Methods.

The multistate Bennett acceptance ratio (MBAR) method is a prevalent approach for computing the free energies of thermodynamic states. In this work, we introduce BayesMBAR, a Bayesian generalization of the MBAR method. By integration of configurations sampled from thermodynamic states with a prior distribution, BayesMBAR computes a posterior distribution of free energies. Using the posterior distribution, we derive free energy estimations and compute their associated uncertainties. Notably, when a uniform prior distribution is used, BayesMBAR recovers the MBAR's result but provides more accurate uncertainty estimates. Additionally, when prior knowledge about free energies is available, BayesMBAR can incorporate this information into the estimation procedure by using nonuniform prior distributions. As an example, we show that by incorporating the prior knowledge about the smoothness of free energy surfaces, BayesMBAR provides more accurate estimates than the MBAR method. Given MBAR's widespread use in free energy calculations, we anticipate BayesMBAR to be an essential tool in various applications of free energy calculations.

Alchemical Free Energy and Hamiltonian Replica Exchange Molecular Dynamics to Compute Hydrofluorocarbon Isotherms in Imidazolium-Based Ionic Liquids.

Ionic liquids (ILs) have shown promise for applications that leverage differential gas solubility in an IL solvent, e.g., gas separations. Although most available literature provides Henry's law constants, the ability to efficiently estimate full isotherms is important for engineering design calculations. Molecular simulation can be used as a tool to predict full isotherms of gas in ILs. However, particle insertions or deletions in a charge-dense IL medium and the sluggish conformational dynamics of ILs present two sampling challenges for these systems. We therefore devised a method that uses Hamiltonian replica exchange (HREX) molecular dynamics (MD) combined with alchemical free energy calculations to compute full solubility isotherms of two different hydrofluorocarbons (HFCs) in imidazolium-based IL binary mixtures. This workflow is significantly faster than the Gibbs ensemble Monte Carlo (GEMC) simulations which fail to deal with the slow conformational relaxation caused by the sluggish dynamics of ILs. Multiple free energy estimators, including thermodynamic integration, free energy perturbation, and multistate Bennett acceptance ratio method, provided consistent results. Overall, the simulated Henry's law constant, isotherm curvature, and solubility trends match experimental results reasonably well. We close by calculating the full solubility isotherms of two HFCs in IL mixtures that have not been reported in the literature, demonstrating the potential of this method to be used for solubility prediction and setting the stage for future computational screening studies that search for the "best" IL to separate azeotropic HFC mixtures.

Thermodynamic and structural description of relative solubility of the flavonoid rutin by DFT calculations and molecular dynamics simulations

Understanding the solubility of flavonoids is of relevance in view of the broad application such as antioxidant and antitumor agents. The solubility of flavonoid rutin has been experimentally investigated in a homologous series of alcohols, acetone, n-heptane and water solvents. In this work theoretical evaluation of Gibbs free energy of solvation was performed at the Density Functional Theory (DFT) and Molecular Dynamics (MD) levels, with the ωB97x-D functional and multistate Bennett acceptance ratio (MBAR) method, respectively. The degree of solubility of a given solute among various solvents is dictated by solute-solvent interactions, which may be represented by the solvation energy. Our calculated relative energy values due to solvent effects, using quantum chemical methods (PCM and SMD solvent models) do not follow strictly the observed solubility trend. However, the PCM model including solute-solvent dispersion interaction, solute-solvent repulsion interaction and solute cavitation energy contribution to the total energy (named DIS,REP,CAV) yielded an overall agreement with experimental solubility profile, with exception of ethanol solvent. A similar satisfactory accordance with experiment was found for MBresults, with a discrepancy also found for 1-butanol, which deserve further investigation. Analysis of spatial and radial distribution functions (SDF and RDF) for pure solvents and the analysis of the contributions of vdW and electrostatic terms for the total free energy reveal that a balance between dispersive and electrostatic interactions between the rutin and the solvent is necessary to solubilize the flavonoid, which corroborates with the values of these energy contributions to the solvation-free energy calculated by the MBAR method. The use of explicit solvent molecules in DFT calculations is precluded due to the size, flexibility, and many sites for solvent interactions of rutin molecule. Our results indicate that in general DFT-PCM (DIS,REP,CAV) and MBAR levels of calculations can satisfactorily predict the solvation free energy and solubility trend for flavonoid rutin along a series of solvents, and also provide a reasonable estimate of relative solubilities of other flavonoids, information that can be relevant for further molecular modeling studies.

Coarse-Grained Modeling of Multiple Pathways in Conformational Transitions of Multi-Domain Proteins.

Large-scale conformational transitions in multi-domain proteins are often essential for their functions. To investigate the transitions, it is necessary to explore multiple potential pathways, which involve different intermediate structures. Here, we present a multi-basin (MB) coarse-grained (CG) structure-based Go̅ model for describing transitions in proteins with more than two moving domains. This model is an extension of our dual-basin Go̅ model in which system-dependent parameters are determined systematically using the multistate Bennett acceptance ratio method. In the MB Go̅ model for multi-domain proteins, we assume that intermediate structures may have partial inter-domain native contacts. This approach allows us to search multiple transition pathways that involve distinct intermediate structures using the CG molecular dynamics (MD) simulations. We apply this scheme to an enzyme, adenylate kinase (AdK), which has three major domains and can move along two different pathways. Using the optimized mixing parameters for each pathway, AdK shows frequent transitions between the Open, Closed, and the intermediate basins and samples a wide variety of conformations within each basin. The explored multiple transition pathways could be compared with experimental data and examined in more detail by atomistic MD simulations.

Extension of the Variational Free Energy Profile and Multistate Bennett Acceptance Ratio Methods for High-Dimensional Potential of Mean Force Profile Analysis.

We redevelop the variational free energy profile (vFEP) method using a cardinal B-spline basis to extend the method for analyzing free energy surfaces (FESs) involving three or more reaction coordinates. We also implemented software for evaluating high-dimensional profiles based on the multistate Bennett acceptance ratio (MBAR) method which constructs an unbiased probability density from global reweighting of the observed samples. The MBAR method takes advantage of a fast algorithm for solving the unbinned weighted histogram (UWHAM)/MBAR equations which replaces the solution of simultaneous equations with a nonlinear optimization of a convex function. We make use of cardinal B-splines and multiquadric radial basis functions to obtain smooth, differentiable MBAR profiles in arbitrary high dimensions. The cardinal B-spline vFEP and MBAR methods are compared using three example systems that examine 1D, 2D, and 3D profiles. Both methods are found to be useful and produce nearly indistinguishable results. The vFEP method is found to be 150 times faster than MBAR when applied to periodic 2D profiles, but the MBAR method is 4.5 times faster than vFEP when evaluating unbounded 3D profiles. In agreement with previous comparisons, we find the vFEP method produces superior FESs when the overlap between umbrella window simulations decreases. Finally, the associative reaction mechanism of hammerhead ribozyme is characterized using 3D, 4D, and 6D profiles, and the higher-dimensional profiles are found to have smaller reaction barriers by as much as 1.5 kcal/mol. The methods presented here have been implemented into the FE-ToolKit software package along with new methods for network-wide free energy analysis in drug discovery.

Replica-Exchange Umbrella Sampling Combined with Gaussian Accelerated Molecular Dynamics for Free-Energy Calculation of Biomolecules.

We have developed an enhanced conformational sampling method combining replica-exchange umbrella sampling (REUS) with Gaussian accelerated molecular dynamics (GaMD). REUS enhances the sampling along predefined reaction coordinates, while GaMD accelerates the conformational dynamics by adding a boost potential to the system energy. The method, which we call GaREUS (Gaussian accelerated replica-exchange umbrella sampling), enhances the sampling more efficiently than REUS or GaMD, while the computational resource for GaREUS is the same as that required for REUS. The two-step reweighting procedure using the multistate Bennett acceptance ratio method and the cumulant expansion for the exponential average is applied to the simulation trajectories for obtaining the unbiased free-energy landscapes. We apply GaREUS to the calculations of free-energy landscapes for three different cases: conformational equilibria of N-glycan, folding of chignolin, and conformational change of adenyl kinase. We show that GaREUS speeds up the convergences of free-energy calculations using the same amount of computational resources as REUS. The free-energy landscapes reweighted from the trajectories of GaREUS agree with previously reported ones. GaREUS is applicable to free-energy calculations of various biomolecular dynamics and functions with reasonable computational costs.

Q-Canonical Monte Carlo Sampling for Modeling the Linkage between Charge Regulation and Conformational Equilibria of Peptides.

The overall charge content and the patterning of charged residues have a profound impact on the conformational ensembles adopted by intrinsically disordered proteins. These parameters can be altered by charge regulation, which refers to the effects of post-translational modifications, pH-dependent changes to charge, and conformational fluctuations that modify the pKa values of ionizable residues. Although atomistic simulations have played a prominent role in uncovering the major sequence-ensemble relationships of IDPs, most simulations assume fixed charge states for ionizable residues. This may lead to erroneous estimates for conformational equilibria if they are linked to charge regulation. Here, we report the development of a new method we term q-canonical Monte Carlo sampling for modeling the linkage between charge regulation and conformational equilibria. The method, which is designed to be interoperable with the ABSINTH implicit solvation model, operates as follows: For a protein sequence with n ionizable residues, we start with all 2n charge microstates and use a criterion based on model compound pKa values to prune down to a subset of thermodynamically relevant charge microstates. This subset is then grouped into mesostates, where all microstates that belong to a mesostate have the same net charge. Conformational distributions, drawn from a canonical ensemble, are generated for each of the charge microstates that make up a mesostate using a method we designate as proton walk sampling. This method combines Metropolis Monte Carlo sampling in conformational space with an auxiliary Markov process that enables interconversions between charge microstates along a mesostate. Proton walk sampling helps identify the most likely charge microstate per mesostate. We then use thermodynamic integration aided by the multistate Bennett acceptance ratio method to estimate the free energies for converting between mesostates. These free energies are then combined with the per-microstate weights along each mesostate to estimate standard state free energies and pH-dependent free energies for all thermodynamically relevant charge microstates. The results provide quantitative estimates of the probabilities and preferred conformations associated with every thermodynamically accessible charge microstate. We showcase the application of q-canonical sampling using two model systems. The results establish the soundness of the method and the importance of charge regulation in systems characterized by conformational heterogeneity.

Fast Solver for Large Scale Multistate Bennett Acceptance Ratio Equations.

The multistate Bennett acceptance ratio method (MBAR) and unbinned weighted histogram analysis method (UWHAM) are widely employed approaches to calculate relative free energies of multiple thermodynamic states that gain statistical precision by employing free energy contributions from configurations sampled at each of the simulated λ states. With the increasing availability of high throughput computing resources, a large number of configurations can be sampled from hundreds or even thousands of states. Combining sampled configurations from all states to calculate relative free energies requires the iterative solution of large scale MBAR/UWHAM equations. In the current work, we describe the development of a fast solver to iteratively solve these large scale MBAR/UWHAM equations utilizing our previous findings that the MBAR/UWHAM equations can be derived as a Rao-Blackwell estimator. The solver is implemented and distributed as a Python module called FastMBAR. Our benchmark results show that FastMBAR is more than 2 times faster than the widely used solver pymbar, when it runs on a central processing unit (CPU) and more than 100 times faster than pymbar when it runs on a graphical processing unit (GPU). The significant speedup achieved by FastMBAR running on a GPU is useful not only for solving large scale MBAR/UWHAM equations but also for estimating uncertainty of calculated free energies using bootstrapping where the MBAR/UWHAM equations need to be solved multiple times.

Evaluation of the SAFT-γ Mie force field with solvation free energy calculations

Solvation free energies can be essential in the process of evaluating and developing force fields. In addition, they can be used to obtain a diversity of other thermodynamic properties such as Henry's law constants and partition coefficients, for instance. In this paper, we study the solvation of asphaltene-like molecules (polyaromatic hydrocarbons) in both aqueous and organic solvents by using a coarse-grained model known as the SAFT-γ Mie force field. This model relies on a top-down approach in which the force-field parameters are obtained using an equation of state. The use of solvation free energy calculations to evaluate this force field can help in improving the model and increasing the scale in which these simulations can be applied. The results presented here were obtained by carrying out molecular dynamics simulations, with the expanded ensemble method being applied to sample an alchemical path. For this, we employed a softcore variant of the Mie potential. From the output of these simulations, we estimated free-energy differences by means of the Multistate Bennett Acceptance Ratio method. The results with organic solvents exhibited low absolute deviations from experimental data. In turn, hydration free energy calculations required a binary interaction parameter to be estimated from molecular dynamics data, which displays a flaw in the top-down parameterization approach. Fortunately, however, a single binary parameter could be used for all pairs of polyaromatic hydrocarbons with water, proving that the SAFT-γ Mie force field exhibits a suitable transferability property that deserves further investigation.

A comparative study of methods to compute the free energy of an ordered assembly by molecular simulation

We present a comparative study of methods to compute the absolute free energy of a crystalline assembly of hard particles by molecular simulation. We consider all combinations of three choices defining the methodology: (1) the reference system: Einstein crystal (EC), interacting harmonic (IH), or r(-12) soft spheres (SS); (2) the integration path: Frenkel-Ladd (FL) or penetrable ramp (PR); and (3) the free-energy method: overlap-sampling free-energy perturbation (OS) or thermodynamic integration (TI). We apply the methods to FCC hard spheres at the melting state. The study shows that, in the best cases, OS and TI are roughly equivalent in efficiency, with a slight advantage to TI. We also examine the multistate Bennett acceptance ratio method, and find that it offers no advantage for this particular application. The PR path shows advantage in general over FL, providing results of the same precision with 2-9 times less computation, depending on the choice of a common reference. The best combination for the FL path is TI+EC, which is how the FL method is usually implemented. For the PR path, the SS system (with either TI or OS) proves to be most effective; it gives equivalent precision to TI+FL+EC with about 6 times less computation (or 12 times less, if discounting the computational effort required to establish the SS reference free energy). Both the SS and IH references show great advantage in capturing finite-size effects, providing a variation in free-energy difference with system size that is about 10 times less than EC. This result further confirms previous work for soft-particle crystals, and suggests that free-energy calculations for a structured assembly be performed using a hybrid method, in which the finite-system free-energy difference is added to the extrapolated (1/N→0) absolute free energy of the reference system, to obtain a result that is nearly independent of system size.

Minimum Free Energy Path of Ligand-Induced Transition in Adenylate Kinase

Large-scale conformational changes in proteins involve barrier-crossing transitions on the complex free energy surfaces of high-dimensional space. Such rare events cannot be efficiently captured by conventional molecular dynamics simulations. Here we show that, by combining the on-the-fly string method and the multi-state Bennett acceptance ratio (MBAR) method, the free energy profile of a conformational transition pathway in Escherichia coli adenylate kinase can be characterized in a high-dimensional space. The minimum free energy paths of the conformational transitions in adenylate kinase were explored by the on-the-fly string method in 20-dimensional space spanned by the 20 largest-amplitude principal modes, and the free energy and various kinds of average physical quantities along the pathways were successfully evaluated by the MBAR method. The influence of ligand binding on the pathways was characterized in terms of rigid-body motions of the lid-shaped ATP-binding domain (LID) and the AMP-binding (AMPbd) domains. It was found that the LID domain was able to partially close without the ligand, while the closure of the AMPbd domain required the ligand binding. The transition state ensemble of the ligand bound form was identified as those structures characterized by highly specific binding of the ligand to the AMPbd domain, and was validated by unrestrained MD simulations. It was also found that complete closure of the LID domain required the dehydration of solvents around the P-loop. These findings suggest that the interplay of the two different types of domain motion is an essential feature in the conformational transition of the enzyme.

Conformational Transition Pathways of Adenylate Kinase Explored by the String Method

Large-scale conformational changes in proteins involve barrier-crossing transitions on the complex free energy surfaces of high-dimensional space. Here we show that, by combining the on-the-fly string method and the multi-state Bennett acceptance ratio (MBAR) method, the free energy profile of a conformational transition pathway in Escherichia coli adenylate kinase can be accurately determined. The minimum free energy paths (MFEPs) of adenylate kinase were searched by the on-the-fly string method in 20-dimensional space spanned by the 20 largest-amplitude principal modes, and the free energy and various kinds of average physical quantities along the MFEPs were successfully evaluated by the MBAR method. The influence of ligand binding on the MFEPs was characterized in terms of rigid-body motions of the LID and AMPbd domains. It was found that the LID domain was able to partially close without the ligand, while the closure of the AMPbd domain required the ligand binding. The transition state ensemble of the ligand bound form was identified as those structures characterized by highly specific binding of the ligand to the AMPbd domain, and was validated by unrestrained MD simulations. It was also found that complete closure of the LID domain required the dehydration of solvents around the P-loop. These findings suggest that the interplay of the two different types of domain motion is an feature in the conformational transition of the enzyme.

Using multistate free energy techniques to improve the efficiency of replica exchange accelerated molecular dynamics.

Replica exchange accelerated molecular dynamics (REXAMD) is a method that enhances conformational sampling while retaining at least one replica on the original potential, thus avoiding the statistical problems of exponential reweighting. In this article, we study three methods that can combine the data from the accelerated replicas to enhance the estimate of properties on the original potential: weighted histogram analysis method (WHAM), pairwise multistate Bennett acceptance ratio (PBAR), and multistate Bennett acceptance ratio (MBAR). We show that the method that makes the most efficient use of equilibrium data from REXAMD simulations is the MBAR method. This observation holds for both alchemical free energy and structural observable prediction. The combination of REXAMD and MBAR should allow for more efficient scaling of the REXAMD method to larger biopolymer systems.